Kernel Density Estimation

Installing maptools

install.packages("maptools", repos = "https://packagemanager.posit.co/cran/2023-10-13")

pacman::p_load(maptools, sf, raster, spatstat, tmap, tidyverse)

Spatial Data Wrangling

Importing the spatial data

In this section, st_read() of sf package will be used to import these three geospatial data sets into R.

childcare_sf <- st_read("data/ChildCareServices.geojson") %>%
  st_transform(crs = 3414)

Reading layer `ChildCareServices' from data source 
  `D:\tskam\IS415-GAA\In-class_Ex\In-class_Ex03\data\ChildCareServices.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

#sg_sf <- st_read(dsn = "data", layer="CostalOutline")
mpsz_sf <- st_read(dsn = "data", layer="MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `D:\tskam\IS415-GAA\In-class_Ex\In-class_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

plot(mpsz_sf)

Creating coastal outline

sg_sf <- mpsz_sf %>%
  st_union()

plot(sg_sf)

Notice that except childcare_sf, both mpsz_sf and sg_sf do not have proper crs information.

DIY: Using the method you learned in Lesson 2, assign the correct crs to mpsz_sf and sg_sf simple feature data frames.

DIY: If necessary, changing the referencing system to Singapore national projected coordinate system.

Mapping the geospatial data sets

After checking the referencing system of each geospatial data data frame, it is also useful for us to plot a map to show their spatial patterns.

DIY: Using the mapping methods you learned in Hands-on Exercise 3, prepare a map as shown below.

Notice that all the geospatial layers are within the same map extend. This shows that their referencing system and coordinate values are referred to similar spatial context. This is very important in any geospatial analysis.

Alternatively, we can also prepare a pin map by using the code chunk below.

tmap_mode('view')
tm_shape(childcare_sf)+
  tm_dots()

tmap_mode('plot')

Notice that at the interactive mode, tmap is using leaflet for R API. The advantage of this interactive pin map is it allows us to navigate and zoom around the map freely. We can also query the information of each simple feature (i.e. the point) by clicking of them. Last but not least, you can also change the background of the internet map layer. Currently, three internet map layers are provided. They are: ESRI.WorldGrayCanvas, OpenStreetMap, and ESRI.WorldTopoMap. The default is ESRI.WorldGrayCanvas.

Reminder: Always remember to switch back to plot mode after the interactive map. This is because, each interactive mode will consume a connection. You should also avoid displaying ecessive numbers of interactive maps (i.e. not more than 10) in one RMarkdown document when publish on Netlify.

Geospatial Data wrangling

Creating ppp objects: sf method

childcare_ppp <- as.ppp(childcare_sf)

summary(childcare_ppp)

Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1925 character character 

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units

plot(childcare_ppp)

Notice the warning message about duplicates. In spatial point patterns analysis an issue of significant is the presence of duplicates. The statistical methodology used for spatial point patterns processes is based largely on the assumption that process are simple, that is, that the points cannot be coincident.

Handling duplicated points

We can check the duplication in a ppp object by using the code chunk below.

any(duplicated(childcare_ppp))

[1] FALSE

To count the number of co-indicence point, we will use the multiplicity() function as shown in the code chunk below.

multiplicity(childcare_ppp)

If we want to know how many locations have more than one point event, we can use the code chunk below.

sum(multiplicity(childcare_ppp) > 1)

[1] 0

The output shows that there are 128 duplicated point events.

To view the locations of these duplicate point events, we will plot childcare data by using the code chunk below.

tmap_mode('view')
tm_shape(childcare_sf) +
  tm_dots(alpha=0.4, 
          size=0.05)

tmap_mode('plot')

Challenge: Do you know how to spot the duplicate points from the map shown above?

There are three ways to overcome this problem. The easiest way is to delete the duplicates. But, that will also mean that some useful point events will be lost.

The second solution is use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

The third solution is to make each point “unique” and then attach the duplicates of the points to the patterns as marks, as attributes of the points. Then you would need analytical techniques that take into account these marks.

The code chunk below implements the jittering approach.

childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

DIY: Using the method you learned in previous section, check if any dusplicated point in this geospatial data.

any(duplicated(childcare_ppp_jit))

[1] FALSE

Creating owin object: sf method

sg_owin <- as.owin(sg_sf)

The ouput object can be displayed by using plot() function

plot(sg_owin)

and summary() function of Base R.

summary(sg_owin)

Combining point events object and owin object

In this last step of geospatial data wrangling, we will extract childcare events that are located within Singapore by using the code chunk below.

childcareSG_ppp = childcare_ppp[sg_owin]

The output object combined both the point and polygon feature in one ppp object class as shown below.

summary(childcareSG_ppp)

DIY: Using the method you learned in previous exercise, plot the newly derived childcareSG_ppp as shown below.

First-order Spatial Point Patterns Analysis

In this section, you will learn how to perform first-order SPPA by using spatstat package. The hands-on exercise will focus on:

• deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes,
• performing Confirmatory Spatial Point Patterns Analysis by using Nearest Neighbour statistics.

Kernel Density Estimation

In this section, you will learn how to compute the kernel density estimation (KDE) of childcare services in Singapore.

Computing kernel density estimation using automatic bandwidth selection method

The code chunk below computes a kernel density by using the following configurations of density() of spatstat:

• bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl().
  
• The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”.
  
• The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.

kde_childcareSG_bw <- density(childcareSG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

The plot() function of Base R is then used to display the kernel density derived.

plot(kde_childcareSG_bw)

The density values of the output range from 0 to 0.000035 which is way too small to comprehend. This is because the default unit of measurement of svy21 is in meter. As a result, the density values computed is in “number of points per square meter”.

Before we move on to next section, it is good to know that you can retrieve the bandwidth used to compute the kde layer by using the code chunk below.

bw <- bw.diggle(childcareSG_ppp)
bw

   sigma 
295.4419 

Rescalling KDE values

In the code chunk below, rescale() is used to covert the unit of measurement from meter to kilometer.

childcareSG_ppp.km <- rescale(childcareSG_ppp, 
                              1000, 
                              "km")

Now, we can re-run density() using the resale data set and plot the output kde map.

kde_childcareSG.bw <- density(childcareSG_ppp.km,
                              sigma=bw.diggle,
                              edge=TRUE,
                              kernel="gaussian")
plot(kde_childcareSG.bw)

Notice that output image looks identical to the earlier version, the only changes in the data values (refer to the legend).

Working with different automatic badwidth methods

Beside bw.diggle(), there are three other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().

Let us take a look at the bandwidth return by these automatic bandwidth calculation methods by using the code chunk below.

 bw.CvL(childcareSG_ppp.km)

  sigma 
4.54311 

bw.scott(childcareSG_ppp.km)

 sigma.x  sigma.y 
2.159749 1.396455 

bw.ppl(childcareSG_ppp.km)

    sigma 
0.3897017 

bw.diggle(childcareSG_ppp.km)

    sigma 
0.2954419 

Baddeley et. (2016) suggested the use of the bw.ppl() algorithm because in ther experience it tends to produce the more appropriate values when the pattern consists predominantly of tight clusters. But they also insist that if the purpose of once study is to detect a single tight cluster in the midst of random noise then the bw.diggle() method seems to work best.

The code chunk beow will be used to compare the output of using bw.diggle and bw.ppl methods.

kde_childcareSG.ppl <- density(childcareSG_ppp.km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

Working with different kernel methods

By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics.

The code chunk below will be used to compute three more kernel density estimations by using these three kernel function.

par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")

Fixed and Adaptive KDE

Computing KDE by using fixed bandwidth

Next, you will compute a KDE layer by defining a bandwidth of 600 meter. Notice that in the code chunk below, the sigma value used is 0.6. This is because the unit of measurement of childcareSG_ppp.km object is in kilometer, hence the 600m is 0.6km.

kde_childcareSG_600 <- density(childcareSG_ppp.km,
                               sigma=0.6, 
                               edge=TRUE,
                               kernel="gaussian")
plot(kde_childcareSG_600)

Computing KDE by using adaptive bandwidth

Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. One way to overcome this problem is by using adaptive bandwidth instead.

In this section, you will learn how to derive adaptive kernel density estimation by using density.adaptive() of spatstat.

kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

We can compare the fixed and adaptive kernel density estimation outputs by using the code chunk below.

par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

Converting KDE output into grid object.

The result is the same, we just convert it so that it is suitable for mapping purposes

gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
spplot(gridded_kde_childcareSG_bw)

Converting gridded output into raster

Next, we will convert the gridded kernal density objects into RasterLayer object by using raster() of raster package.

kde_childcareSG_bw_raster <- raster(gridded_kde_childcareSG_bw)

Let us take a look at the properties of kde_childcareSG_bw_raster RasterLayer.

kde_childcareSG_bw_raster

class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.419757, 0.2695907  (x, y)
extent     : 2.667538, 56.39644, 15.74872, 50.25633  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -1.293897e-14, 37.27443  (min, max)

Notice that the crs property is NA.

Assigning projection systems

The code chunk below will be used to include the CRS information on kde_childcareSG_bw_raster RasterLayer.

projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster

class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.419757, 0.2695907  (x, y)
extent     : 2.667538, 56.39644, 15.74872, 50.25633  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -1.293897e-14, 37.27443  (min, max)

Notice that the crs property is completed.

Visualising the output in tmap

Finally, we will display the raster in cartographic quality map using tmap package.

tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

Notice that the raster values are encoded explicitly onto the raster pixel using the values in “v”” field.

Comparing Spatial Point Patterns using KDE

In this section, you will learn how to compare KDE of childcare at Ponggol, Tampines, Chua Chu Kang and Jurong West planning areas.

Extracting study area

The code chunk below will be used to extract the target planning areas.

pg = mpsz_sf[mpsz_sf@data$PLN_AREA_N == "PUNGGOL",]
tm = mpsz_sf[mpsz_sf@data$PLN_AREA_N == "TAMPINES",]
ck = mpsz_sf[mpsz_sf@data$PLN_AREA_N == "CHOA CHU KANG",]
jw = mpsz_sf[mpsz_sf@data$PLN_AREA_N == "JURONG WEST",]

pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")

Plotting target planning areas

par(mfrow=c(2,2))
plot(pg, main = "Ponggol")

plot(tm, main = "Tampines")

plot(ck, main = "Choa Chu Kang")

plot(jw, main = "Jurong West")